



2. Government Accession No. 


3. Recipient’s Catalog No. 


1. Report No. 

NASA TR R-454 

4, Title and Subtitle I ^ 


PRIMER VECTOR THEORY AND APPLICATIONS 


7. Author(s) 8. 


Report Date 

November 1975 

Performing Organization Code 

JSC-09259 


Performing Organization Report No. 


Donald J. Jezewski 


9. Performing Organization Name and Address 

Lyndon B. Johnson Space Center 
Houston, Texas 77058 


1 2. Sponsoring Agency Name and Address 


JSC S-439 

10. Work Unit No. 

986-16-20-00-72 

11. Contract or Grant No. 


13. Type of Report and Period Covered 

Technical Report 


National Aeronautics and Space Administration 
Washington, D. C. 20546 


15. Supplementary Notes 


16. Abstract 


Sponsoring Agency Code 


A method developed to compute two- body, optimal, N- impulse trajectories is presented. The 
necessary conditions established define the gradient structure of the primer vector and its 
derivative for any set of boundary conditions and any number of impulses. Inequality constraints, 
a conjugate gradient iterator technique, and the use of a penalty function also are discussed. 


17. Key Words (Suggested by Author(s)} 

‘Trajectory Optimization 
‘ Pontryagin Principle 
‘Calculus of Variations 
‘ Lagrange Multipliers 


* Transfer Orbits 
' Impulses 


19. Security Classif. (of this report) 

Unclassified 


18. Distribution Statement 


STAR Subject Category: 
13 (Astrodynamics) 


20. Security Classif. (of this page) 

Unclassified 


21. No. of Pages 

47 


22. Price* 

$3.75 


*For sale by the National Technical Information Service, Springfield. Virginia 22 161 


I I 1 . £ £ 


i I i I M I i I 


1 1 1 1 1 










[ 


1 I i I < i I I I M I I I i I 1 I I 



CONTENTS 


Section 

SUMMARY 

INTRODUCTION 

DEFINITION 

LAWDEN'S NECESSARY CONDITIONS 

CRITERIA FOR AN ADDITIONAL IMPULSE . . . 

Adjoint Equation 

State Transition Matrix 

COST FUNCTION 

Fixed End Conditions and Fixed Transfer Time 

General Case 

EXAMPLE PROBLEMS 

Fixed-End- Condition Transfer 

Orbit- to- Orbit Transfer 

Generalized Rendezvous 

INEQUALITY CONSTRAINTS IN PRIMER OPTIMAL 
N-IMPULSE SOLUTIONS 

Periapsis Constraint 

Lighting Constraint 

CONCLUDING REMARKS 

REFERENCES 


Page 

1 

1 

2 

3 

11 

11 

13 

14 
14 
24 
28 
28 
29 
31 

34 

34 

41 

42 

43 


111 


1 I i ;; 


I I I I M I i E I s 


I I i 



FIGURES 


Figure 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 


Thrust-vector orientation 

Bang-bang control 

Switch function 

(a) Finite thrust 

(b) Impulsive thrust 

Switch function, impulsive thrust, t = tj 

Primer magnitude 

Reference and perturbed trajectories (fixed conditions) , 
where fixed transfer time ^ ' ^2 ~ H 

Primer magnitude characteristics 

Primer magnitude (nonoptimal) characteristics .... 

Reference and perturbed trajectories (general case) . 

Transfer between two orbits 

Optimal three-impulse rendezvous 

Radius-of-periapsis inequality constraint 

Single-wall penalty function 

Double-wall penalty function (cost well) 

Coordinate system for lighting inequality constraint 


Page 

3 

3 


9 

9 

9 

10 

15 

19 

23 

25 

30 

33 

35 

36 
36 
42 


f 


[ 


r 


iv 


I ! I I i . 





I 



PRIMER VECTOR THEORY AND APPLICATIONS 

By Donald J . Jezewski 
Lyndon B. Johnson Space Center 

SUMMARY 


The results of earlier studies establishing the criteria for an optimal impulsive 
trajectory and for improving a reference impulsive trajectory are used in this re- 
port to develop the necessary conditions for a general, two-body, optimal, N-impulse 
solution. A general differential cost function is developed which defines the gradient 
structure for any set of boundary conditions when the applicable constraints are 
specified . 

Example problems are presented to illustrate the generality of the solution; 
these problems concern fixed-end-condition transfer , orbit-to-orbit transfer , and 
generalized rendezvous . 

Use of a penalty function approach (cost well) enables the establishment of 
inequality constraints on both state and control variables . The periapsis constraint 
and a lighting constraint are discussed specifically . 


INTRODUCTION 


In the mid to late 1960's, the solution of an optimal trajectory based on replacing 
the thrust with an impulse received considerable attention as a potentially rapid and 
reliable method for performing mission design studies. This renewed interest in 
impulse solutions was prompted by the extensive work of Lawden (ref. 1) . Lawden 
developed the necessary conditions for an optimal impulsive trajectory by examining 
the limiting conditions on an optimal finite thrust solution wherein the thrust magni- 
tude was unconstrained but bounded between a maximum and a minimum value. The 
results he obtained are known today as Lawden' s necessary conditions for an optimal 
impulse trajectory . These results specify the conditions that must be satisfied by the 
primer vector and its derivative on an impulsive trajectory that is considered optimal . 

In early 1968, Lion and Handelsman (ref. 2) established the criteria whereby a 
reference impulsive trajectory could be improved (i.e. , whereby a decrease in the 
cost function or the sum of the magnitudes of the applied impulses could be accom- 
plished). These criteria established the conditions under which an additional inter- 
mediate impulse or an initial or final coast would improve the solution. Using Lawden's 
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necessary conditions , they developed a first-order differential (a gradient vector) 
of the cost function with respect to the intermediate position vectors and time . With 
this development , the mechanism for systematically determining the minimum cost 
for an N-impulse trajectory was created. 

In late 1968, Jezewski and Rozendaal perceived the opportunity to generate an 
efficient method for computing two-body, optimal N-impulse trajectories by combin- 
ing Lion and Handelsman’s gradient vector with a conjugate gradient iterator. A 
computer program known as the Optimal, Multi, Impulse Rendezvous (OMIR) program 
was produced. Fixed-time , N-impulse (as many as six impulses), optimal trajectories 
(ref. 3) could be efficiently and quickly generated with this algorithm. 

Boundary conditions (open time, initial and terminal constraints, rendezvous, 
etc . ) were extended with the generation of a general differential cost function . This 
function defines the gradient structure of the cost function for any set of boundary 
conditions when the applicable constraints are specified. Inequality constraints 
(initial, intermediate, and terminal) on both state and control variables were imple 
mented into the OMIR program by a penalty function approach (cost well) (ref. 4) . 
By this means, completely general, two-body, N-impulse, optimal trajectories can be 
generated for any set of constraints that can be expressed mathematically . 

This report consists of material used by the author in a series of lectures on 
primer vector theory at the NASA Lyndon B . Johnson Space Center in the fall of 1970 
and the spring of 1971 . These lectures were designed to acquaint the users of the 
OMIR program with the necessary background material on primer vector theory and 
its application to trajectory problems. Two approaches are used to present the 
material: the classical calculus of variations and Pontryagin's maximum principle. 

The approach used for each topic was determined by the ease and clarity with which 
the material for that topic could be presented . 

DEFINITION 


In a rectangular Cartesian coordinate system, the primer vector P is defined 
as a vector that is parallel to the thrust vector T (fig. 1). The symbol P does not 

represent the thrust direction cosine vector B , although the optimum value of this 
vector 8 can be defined in terms of the primer vector as 8^^^ = P/|P| . For an 

optimum transfer problem , the solution requires determining the thrust (magnitude 
and direction) , as a function of time , that minimizes some performance index while 
satisfying a boundary condition (BC) set. Although the thrust magnitude may not be 
continuous throughout the solution, the magnitude will usually be bounded (fig. 2) . 

Other names associated with the primer vector and its derivative are Lagrange 
multipliers, costate vectors, adjoint variables, and influence coefficients. To better 
understand the primer vector and its derivative, their origin, and their use in de- 
termining optimum control, it will be helpful to examine an optimal trajectory for a 
mass particle moving in a general gravitational field , in a vacuum , for which the 
thrust-vector time history is not predetermined. If the results obtained from this 
solution are examined for a thrust that is impulsive ( ] T j “ or the ratio of the burn 
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time to the solution time 0) , the necessary conditions on the primer P and 

its derivative P will be obtained for an optimal impulsive trajectory. 


LAWDEN'S NECESSARY CONDITIONS 


The following material is an interpretation of reference 1 , supplemented by 
references 5 to 7. The state equations for the trajectory problem are 


m 

(1) 

> 

II 

(2) 

m = -p 

(3) 


wherein R, V , and G are the position, velocity, and gravity vectors, respectively . 
The thrust magnitude is specified by the product pc , where p is the mass flow rate , 
c is the exhaust velocity of the engine , and m is the mass of the vehicle . The 
direction cosine vector C is constrained by the relationship 


T 

C £ = 1 


(4) 
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The mass flow rate is constrained by the relationship 


P(Pmax-P) 


( 5 ) 


wherein a is known as a slack variable. The initial conditions (R(0),V(0) ,m(0)) 
and the control variables over the solution (8 (t) ,p (t) ,a(t)) are needed to obtain a 
solution of this system of equations (R(t) ,V(t) ,m(t)) , as shown in the following 
diagram . 



I 


It is desirable to minimize a performance index J 




( 6 ) 


subject to the differential constraints (eqs . (1) to (3)) and the algebraic constraints 
(eqs. (4) and (5)). The subscript F refers to an evaluation at the final time. Thus, 
the Hamiltonian function is |T 

H = , g] + qT'v - nP 


wherein the vectors P and Q and the scalars q , , and ^2 sre time varying 

Lagrange multipliers . 
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The necessary conditions for an optimal trajectory are as follows. 





0 = 


0 = 


0 = 


= -o^ 
av ^ 


dH 

'dR 


T ^ 
3R 


-I? = P4p'^£ 


3m 


m 


dH 

32 


3H 

ap 


aa 


= + 2n 

m ^1 

= + T] + p„/p 

m ' ^2 ^^max 

= -2^2“ 



( 8 ) 

(9) 

( 10 ) 

( 11 ) 

( 12 ) 

(13) 


From equation (13) , ^2 ^ ® ^ ® ^2 " “ = If 


a = 0, 


P = 0 

or 

p - p 


max 


(14) 


(Seeeq. (5).) If a#0 and P 2 = II’ ° P ^ ^max' Thus , the solution may have 
the following three types of arcs . 

1. Maximum thrust (MT) 

2. Null thrust (NT), free-fall arc 

3. Intermediate thrust (IT), singular arc 
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From equation (11) , the following conditions are apparent. 

1. If p ^ 0 and P and E are parallel. 

2. If p = 0 and = 0, E is indeterminate. ^ 

3 . If = 0 and p 0 , P vanishes . 

Thus, P and E are parallel; however , are they parallel in the same or in opposing 

directions? The Weierstrass condition is used to resolve this problem . 

The Weierstrass E function is defined in reference 5 as 

I 

9 F ^ 

E = F(X,X,t) - F(X*,X*,t) - ^(X*.X*,t)(X - X*) > 0 (15) 


wherein 


F = T^(F - X) 



P 


V 

T = 

Q 

X = 

R 

1 

n 

1 1 

m 


(16) 


f 


and wherein * indicates optimal trajectory and F represents the right-hand side of 
the state equations. Because the state equations are satisfied on the optimal as well 
as on the nearby trajectory, F(X,X,t) = F(X*,X*,t) = 0 and equation (15) becomes 


E = T^(X - X*) > 0 


L 

(17) 


Substituting for ^ and X and canceling like terms results in 




(18) 


Now , examine this equation for the three types of arcs . 
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1 . For the maximum thrust arc (p = P ) , if the optimum mass flow rate 
P* = Pj^g^> equation (18) reduces to 


T 

PC 


T 

> P C* 


( 19 ) 


T 

This condition will be satisfied for all values of C* if and only if P C takes on its 
maximum value for all variations of C . Hence, P must be alined with C or C = P/p 
where p = | P | and 


Ti<^|P| (20) 

c T 

2. For the null thrust arc (p = 0) , equation (18) reduces to q ^ - P C*. 

For a variable C*, the maximum value of the right-hand side is obtained when P 

T 

and C* are alined, or (P C*) =p. Thus, 

max ^ 


n > 


cp 

m 


( 21 ) 


3. For the intermediate thrust arc (p > p* or p < p*), if C = C*, the in- 
equality in equation (18) is not possible for all permissible values of p*. However, 
equality is possible for 


q = (22) 

With this value for q , a further requirement is q > -P^fi* or n > -p. But eoua- 

^ m * — ^ 

tion (22) implies that q < -p and that equality exists only when P and C are 
alined. Thus, for IT arcs, P and C are alined and 


H = 



(23) 
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The Weierstrass condition is summarized as follows. 

1 , If p / 0 and p 0 , thrust is in the direction of the primer . 

2 . If $ is defined as 



(24) 


then it is necessary that $ > 0 for an MT arc, $ < 0 for an NT arc, and $ = 0 for 
an IT arc . 

Consider now corners separating two types of arcs . The parameters P , Q , 
T) , and H are known to be continuous . If H is to be continuous , then 

^p - riP = p$ must be continuous. When p is discontinuous, or the type of arc 

changes the switch function, $ e 0. Taking the time derivative of $, one 

obtains $ = ^ + - p . Substituting from equation (10) for q , one obtains 

m z 
m 


$ - 


cp 

m 


W 


(25) 


On an NT arc, m is constant. Integrating equation (25) , one obtains 


$ = ^ + constant 
m 


(26) 


y 

L 


Consider the computation of equation (12) , ^2 “ g 28 

" ' max ' 

MT arc or an NT arc , and P 2 “ ® 

Concerning impulsive thrust, if p ~ or t^/t,p 0, the maximum thrust is 

replaced by an impulse of negligible time duration . Under these conditions , R , G , 
and 9G/9R are continuous; V is discontinuous at impulse times. Hence, P, Q, 
and Q are continuous; Q is discontinuous (involves velocities) . 

Examine the switch function $ shown in figures 3(a) and 3(b). At the time of 
application of an impulse tj, $ = $ = 0 and $(tj) = maximum. Another possibility 


f 
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(a) Finite thrust . (b) Impulsive thrust . 

Figure 3.- Switch function. 



Figure 4.- Switch function, impulsive 
thrust, t = tj- 


is shown in figure 4 . To prove that $ = $ = 0 at t = t^ for an optimal impulse , con 
sider the time derivative of the Hamiltonian function . 


dH _ an _T3G 

df- P W (27) 
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because the right-hand side of equation (27) is continuous and remains finite. Hence, 
from equation (28) , the Hamiltonian function is continuous across the impulse. But, 
on an optimal trajectory, from equations (7) , (8) , and (24) 

r 

H = P^G - p'^V + p$ (29) 


But p$ = 0 on an NT arc or on an IT arc; hence, p$ must vanish on both sides of 

T . 

the impulse. Also, because P G is continuous, it follows that 

pTy(t,-) = P^v(t;) (30) ! 


or 

P^(v'^ - V~) = 0 (31) 


Because the directions of the primer and the impulse are known to be alined , or 
- V = yP (proportionality factor y ^ 0) , equation (31) is an indication that 

P^P = 0 (32) 

r 

I 

or that P Ip or p = 0. And, from equation (25), $(t^) = 0. 

Consider an NT arc joining two impulses, as shown in figure 5. If p = a at 
one impulse, then from equation (26) , $ = ^(p - a) . Because $ = 0 at the second 


I 



r 


Figure 5 . - Primer magnitude . 
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impulse, it follows that p = a also at the second impulse. Thus, on an optimal tra- 
jectory, p = a at all impulse points. Between impulses, or on NT arcs, p _< a. 

If p > a, an impulse would be implied and the trajectory would not be optimal. 

To summarize: The following are Lawden's necessary conditions for an optimal 
impulsive trajectory; each must be satisfied. 

1 . P and P are continuous . 

2. p = a at impulse times . 

3. p < a on NT arcs separating impulses. 

4. p = 0 at all interior impulses. 

In these conditions, a is the magnitude of P at one impulse. 

CRITERIA FOR AN ADDITIONAL IMPULSE 


The criteria for an additional impulse are discussed in the following subsections . 


Adjoint Equation 

The material in this subsection is based on reference 8 . Particle motion in a 
central force field is described by the equations 


R = V (33a) 

V = VU(R,t) (33b) 


wherein U is the potential and V is the gradient operator . The variation between 
a perturbed and a reference trajectory is represented as 6; for example, 8R = R' - R. 
Taking variations of equations (33a) and (33b) results in 



(34) 
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wherein I is the identity matrix and G is the gravity gradient matrix; 6 is of the 
form 


a^U(R,t) 

9R. 9R. 

1 ] 


(35) 


r 


wherein i = 1,2,3, and j = 1,2,3. By definition, the system of equations adjoint to 
equation (34) is 


9 

. T 

9 

X 

= -A 

X 


[ 

(36) 


wherein cp and X are adjoint variables. Writing equations (34) and (36) in second- 
order form results in 


6R = G 6R 
X = GX 


(37) 


r 


T T 

Premultiplying the first equation by X and the second equation by 6R and then , . 

subtracting, one obtains ' 

L 


X^6R - 6R^X = x’^G 6R 


6r’^GX 


(38) 


• T 

If X 6V is added and subtracted, then 


rp . . rp. rp .. ^ rp 

X 6V + X 6V - 5R X - X 8V = 0 


(39) 


which may be expressed as the total time derivative 


-( 

dt ' 


T * T 1 

X 6V - X 6R) = 0 


(40) 
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Integrating this equation , one obtains 


A.^6V - = constant (41) 


Equation (41) is known as the adjoint equation. The relationship applies for the in- 
terval between any two impulses and will be used extensively in the analysis . 


State Transition Matrix 

The state variable perturbations at any two times between impulses are related 
by the state transition matrix . 


6R(t) 


6R(t) 


= 9(t,t) 


6V(t) 


6V(x) 


(42) 


wherein 4>(t,x) is the six-by-six matrix of partial derivatives. 




9 


11 


9 


21 


(t.t) 


(t,t) 

922(t’^) 


(43) 


For example, cpj^(t,t) - 9R(t)/8R(x). Note from equation (37) that the adjoint 

variables are transformed in the same manner as perturbations in the position vector; 
therefore , the following equation can be written . 


"X(t)* 


■X(T)‘ 


= 9 (t.t) 


-X(t)_ 


>(t). 


(44) 


At the impulse times t - t^ and t = t 2 , the maximum value of the primer magnitude 
(P = A.) is chosen to be unity; that is , 

M'l) = ^ 

/ \ /i I (45) 

^(‘ 2 ) = ^2 = ^''2/h''2lj 
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The BC equations (eq. (45)) together with equation (44) , when t - t2 and t 
uniquely define X(t^) and X(t2) if the matrix cp^2 is nonsingular. Therefore, 

^ ^12 (^’^) 




(* 2 ) 


9 


11 






( 46 ) 


COST FUNCTION 

The impulsive cost function and its gradient vector under specific and general 
conditions are discussed in the following subsections . 


Fixed End Conditions and Fixed Transfer Time 

Define the cost on an impulsive trajectory J in the following manner (refs. 2 
and 3) . 


J = 



i=l 


AV. 

1 


(47) 


r 


Consider the following reference and perturbed trajectories from an orbit A to an 

orbit B (fig. 6) . At a time t^ ^2^ ’ position vector on the f 

reference trajectory is perturbed by an amount 6R^ . Two Lambert solutions are t 

obtained between the following position vectors and times on the perturbed trajectory . 

1. IBCtp.tjl- [R(t^)+ 6 R^.tJ. 

The velocities on the perturbed trajectory differ from those on the reference trajec- ^ 

toi^ as 

1 . at time t^^ 

2. 6V ~ at time t 

m m 

3 . 6V ^ at time t ^ 

m m 

4. 6V2 at time t2 f 
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Figure 6.- Reference and perturbed tra- 
jectories (fixed conditions) , where 
fixed transfer time T = t 2 ~ t^. 


The cost increment to first order between the perturbed J' trajectory and the 
reference J trajectory is defined as 


dJ = J' - J 


wherein 


J = 





( 48 ) 


( 49 ) 


and 


J’ = 

AV, + 6V, 



1 1 





AV, 



( 50 ) 


Digress temporarily and consider the following problem: To first order, what 
is the difference between the absolute values of a perturbed vector and a reference 
vector? Let df equal this difference in the following equation. 


I 1 I 
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i 


= |X + 6x1 - jxj 


= [^(X + 6X)’^(X + 6X)] - 


= x^x] 


= x^x] 


(X + 6X) (X + 6X) 


®x ^ m') rx ^ w 


|X| 1 + 2 + 


1x1 ) [i + e + 

^ |X| ‘‘x 


T 

6X e. 


wherein C^ = X/|X]. Hence, if 6Rj^ is small. 


AV ^SV ^ 

‘I-’ = |lv j - * ®''m - ®''ra 


Av/sVj 


If equation (45) is used. 


dj = XjTfiVj + 6V_^* - sv^‘ - 
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Evaluating the adjoint equation (eq. (41)) at the boundaries of the two Lambert solu- 
tions of the perturbed trajectory, one obtains 


T • T -T - • -T 

X/6V- - A/6R, = X ^6V - X ^8R 

1111mm mm 

^ i ^6R ^ = :^/6V - X 
mm mm 2222 


(54) 


If Rj^ and R2 are considered fixed , then 


6R^ = 6R2 = 0 


(55) 


Also , X and k are evaluated on the reference trajectory and hence are continuous , 


A = ?. = k 

m m m 


= k = X 
m m m 


(56) 


Using equations (55) and (56) in equation (54) , 


T T 

X- 6V, = X 6V 
11mm 


T + • T + T 

X ^6V - X 6R = K 

mm mm 2 2 


• T 

X 6R 
m m 


(57) 


T T 

Solving the first equation for X 6V. and the second equation for X 6V„, 

XX Za La 

and using the results in equation (53) , one obtains 


■ T + • T 

dJ = X 6R - X 8R + 
mm mm 


6V - 6V 
m m 


- X 


m 


(6V - 6V I 

y m m / 


6V 1 (58) 

m ' 
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Now 


(59) 


and dt has been chosen to be zero. Because the position vector is continuous, 
m 

6R ^ = 6R ~ and equation (58) reduces to 
mm ^ 


6R = dR - R dt 
m m m m 


6R = dR_ - R^ dt 
m m mm 


dJ 





(60) 


r 




This equation is homogeneous in 6V - 8V . If c and the unit vector A are 
defined as 


c 




A = 


6V - 8V 
m m 


(61) 


equation (60) reduces to 


r 


r 



(62) 


The criteria for an additional impulse can now be established, for if > 1. then ^ 

dJ < 0 and J > J' ; or the reference-trajectory cost is greater than the perturbed- 

trajectory cost. The reference-trajectory cost may be improved by applying an 

impulse in the direction of X at the time t . To first order, the greatest decrease 
^ mm 

in cost will be obtained if the intermediate impulse is applied when = l^mlmax 

as in figure 7 . 

f 
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Figure 7.- Primer magnitude 
characteristics. 


If a reference trajectory exists for which > 1 at some time t^ , how much 

should this trajectory be perturbed and in which direction so as to obtain a perturbed 
trajectory with a lower cost? From the transition matrix. 


r 




1 

6R 

m 

= 1 

t ,t\ 
m 11 

6Rj 


6V 

m 


< / 

6V^ 


6R 

m 



SRjl 


1 


^m’^) 



5V 

m 


^ / 


i 


If equation (55) is used in equation (63) , the following relationships are established. 

6V = <(-22<'m-'2>®''2- 


19 
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Eliminating 6Vj and 5 V 2 from this set of equations and solving for 6V^ and 


6V , one obtains 
m 


6V = cpooft 
m ^ 221 m 1 




6R 


m 


(64) 


and 


6V = cp„„ t ,t_ 
m ^22\ m 2 


i)^12 (^m’*2) 


6R 


m 


(65) 


By differencing these two equations and recalling that 6R = 5R = 6R , one 
obtains the following equation . m m m 


6V - 6V 
I m m 




m 


. A 6R 


( 66 ) 


f 


wherein 


A = cp 


22 




(67) 


r 

u. 


The result of solving for SR^ from equation (66) is 


1 ^ 

cn A “1 ni 

6R^ = cA -r — 
m A 

m 


( 68 ) 


As long as c (the magnitude of the intermediate impulse) is sufficiently small to 
ensure that the first-order theory holds , J' < J . 


r 
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Digress temporarily to consider the following problem: What is a better way 
of determining the parameter c? Consider the differential cost function equation 
(eq. (48)). 


dJ 


AV^ + 6Vj 


AVj + AVj - 6V2 



+ c 


(69) 


The result of retaining second -order terms in the expansion of this function is 


dJ 


1 

2IAV 


1 


5Vi^6Vi + 2 


AV, 



2 


" 2 AV^ 



2 6V/AV3 



(70) 


From equations (63) and (68) , 



wherein 



(71) 


(72) 
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By using equation (71) in equation (70) , a quadratic in c can be obtained of the 

2 

form dJ = a^ + a^c + 82 c . Taking the partial derivative of dJ with respect to c 
and solving for c from 9(dJ)/9c = 0 gives 


c 




(73) 


r 


This initial three-impulse trajectory is not necessarily optimum (fig. 8 ). Although 
the primer is continuous (eq. ( 66 )) and is a unit vector at each impulse, a cusp may 
occur in the primer curve because the primer derivative has been neglected . Con- 
sider a differential change in the time t and in the position vector R . The two 

° m m 

Lambert solutions consist of trajectories between the following position vectors and 
times . 


1. [R, ,tj -»■ [R + 6 R , t + dt ] . 

^ 1 ’ 1 ^ m m m m 

2. [R + 6 R_, t + dt ] ^ [R„,t»] . 

m mm m 2 2 

wherein R and t now refer to the initial three-impulse trajectory . The differen- 
m m 

tial cost on this initial three-impulse trajectory is 


dJ = + 


6 V - 6 V 
m m 




(53) 


The adjoint equations evaluated on the two Lambert solutions (considering that 

6 R. = 6 R- = 0, that X is continuous, and that X is not continuous) become 
12m 


U 


x/^ 8 V = X - X ^ 6 R 

11mm mm 


T + • +T + T 

X 6 V - X 6 R = X_ 6 V„ 
mm mm 2 2 


(74) 


Using equation (74) in equation (53) results in 


dJ = X ^ 6 V - X ^ 6 R + 

mm mm 


6 V - 6 V 
m m 


T* J. . +T‘ ^ 

- X 6 V + X 6 R (75) 
mm mm 
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Figure 8.- Primer magnitude (non- 
optimal) characteristics. 


Because is a unit vector along the interior impulse , 


+A. 


m 


T 






(76) 


and the appropriate terms in equation (75) cancel. Removal of these canceled terms 
leaves 


• +T + • -T 

dJ = X ^6R - X 6R 

m m m m 


(77) 


But , to first order , the perturbations 6R and 6R are 

mm 


8R = dR - V ■^dt 
m m m m 


6R = dR - V dt 
m m m m 


(78) 


Using equation (78) in equation (77) results in 


^ - k )dt 

ym m J mym m m myrn 


(79) 
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Because the acceleration is continuous (only a function of position; i. e. , V - -pR/| R 1 . 

T • 

where p is the gravitational parameter) on an impulsive trajectory , may ' 

be added to and subtracted from equation (79) without its value being changed. The 

T T 

Hamiltonian function H on a coasting trajectory is H = P G + Q V, which may be r 

written as 


T • • T 

H = XV - XV 


(80) 


By evaluating equation (80) at t^^ and t^ and using the results in equation (79) , 


the cost function becomes 


m 


dJ = (X ^ - X pdR + 
\ m ml m 


(h ^ - H ]dt 
\ m m / m 


(81) 


Thus, equation (81) illustrates that, on an optimal trajectory (dJ = 0), the primer 
vector derivative X and the Hamiltonian function must be continuous across the 
intermediate impulse . 


General Case 

Consider the derivation of a general differential cost function; apply the ref- 
erence and perturbed trajectories shown in figure 9. The velocities on the perturbed 
trajectory differ from those on the reference trajectory by the following amounts. 


1. 

6Vi 

at the time 



2. 


at the time 



3. 

6V 

at the time 

t 



m 


m 




at the time 



4. 

6V 

t 



m 


m 


5. 


at the time 


V 


+ 

at the time 

+ 


6. 

6V„ 

^0 



The velocities 6V- and 6V„^ are perturbations in the initial and final orbit , re- 
1 ^ 

spectively , that result from the velocity variation V^' - V^ in the initial orbit and 

the velocity variation V„ “ V„' in the final orbit. These variations are caused by 

consideration of perturbations in the initial and final position vectors that in the pre- 
vious analysis were assumed to equal zero (eq. (55)). 
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Figure 9 . - Reference and perturbed 
trajectories (general case). 


The cost increment to first order is dJ = J' - J wherein 


J = 

|AVi| . |AV2 



/ + -\ 


II 

AVj + (6V^ - 6V^ J 

+ 






/ + 

6V 

+ 

AV„ + (6V„ - 6V„ 1 

m 


2(2 2 j 


(82) 


The differential cost function to first order is 


dJ 



8V 


m 





(83) 


wherein, once again, X is chosen to be a unit vector in the direction of the impulse; 
X - AV/|AV1 . Evaluating the adjoint equation (eq. (41)) at the boundaries of two 
conics of the perturbed trajectory (eq. (54)) and using the results in equation (83) , 
one obtains 


• T + T - 
dJ = X^ 6R^ - 6Vj + 


+ - • +T + 

5V - 6V + X 6R 
m m m m 


. -'p _ T' 

X 6R - X 
m m m 


(6V ^ - 6V I - 
I m m / 
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The first two terms deal with perturbations in the initial orbit; the last two terms, 
with perturbations in the final orbit . The primer , because it is computed from the 
reference trajectory, is continuous across the intermediate impulse. Thus, 

X ' = X = X . Because X_ is defined as 

m m m m ^ 


X 


(dV - 6V^ ) 
\ m m / 


m 


5V 


m 


8V 


m 


(85) 


the third and sixth terms in equation (84) cancel. Also, because the position vector p 

at t — t is required to be continuous , the following first-order equations may be 1 ^ 

m 

obtained . 


6R ^ = dR - V ^dt_ 
m m m m 


6R = dR - dt 
m m m m 


(78) 


Therefore, using equations (85) and (78) in equation (84) results in 


dJ 




X ^ idt - 

mm mm/m l L i i 


( 86 ) 


Evaluating the Hamiltonian function (eq . (80)) at t^ and t^ results in 


m 


-T • - • -T 

H = X V - 
m mm mm 


j. +T' . + . +'T> + 

H = X V - X V 
m mm mm 
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By recalling that , equation (87) may be solved for the terms 

* "t'T ^ ~ TP — 

X and X . The resulting expressions may be used in equation (86) to 

obtain 

dJ = x/^6R - X " + (x - X ")^dR 

1 1 1 1 ym m J m 

+ (H^-H'Vt - X„'^8R„ - X„^6V/ (88; 

Vm m/m2222 ^ 

The perturbation in the initial and final position and velocity vectors can be written 
to first order as 


6r/ = dR^ - v/dt^ bR^- = dR^ - V^-dt^ 

bVi- = dV^ - V^-dt^ bv/ = dV^ - v/dt2 


(89) 


By combining equations (89) and (88) and by using the Hamiltonian function evalu- 
ated at t and t = tg, one obtains the general differential cost function. 


dJ 


Hi dti - ^^i^dVi - ^iTdRi) + 

- "2 ■i'2 " - ^^<>''2) * - «n,‘)o'm 


(90) 


Equation (90) can be expressed as 


dJ = VJ dZ 


(91) 


wherein V is the gradient operator and Z is a vector of possible control parameters. 








(92) 
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Equation (90) or (91) defines the general differential cost function. From this ex- 
pression, it can be determined how any two-impulse solution or any two-impulse 
segment of an N-impulse trajectory can be improved by the following changes. 

1. Perturbing the initial variables (tj^,Rj^,Vj) r 

Mr 

2. Perturbing the final variables (t 2 iR 2 >V 2 ) 

3. Perturbing the intermediate variables (t ,R ) 

® mm 

If a particular variable (i.e. , the Zdh) is not allowed to be perturbed and hence is 
fixed, then dZ. = 0 and the general differential cost function is unaffected. 

E 

EXAMPLE PROBLEMS 


Fixed-End-Condition Transfer 

Consider a problem designed to find two-impulse solutions between a fixed 

initial state and a fixed final state for a transfer time that is open or free. Because 

intermediate impulses are excluded , R and t are fixed; therefore, _ 

^ m m ^ 


dR = dt =0 (93) 

m m 


In addition, because the initial position and velocity vector (R^,Vj) and the final 
position and velocity vector (R 2 >V 2 ) are fixed 


dR^ = dR2 = dV^ = dV2 = 0 


(94) 


r 

k 


Using equations (93) and (94) in equation (90) results in a differential cost function 
of 


dJ = dt^ - dt2 


(95) 


Because this problem concerns two-impulse solutions , the Hamiltonian function is 
constant across the trajectory. 


«1 = «2 = " 


(96) 


[ 


28 


I 1 I I 


I I I I I I i I i i I I f 


I 



Finally, the transfer time 6 t is defined as 


6 t = t 2 - t^ (97) 

and its derivative is defined as 

d 6 t = dt 2 - dtj (98) 


Therefore, using equations (98) and (96) in equation (95), the differential cost 
function is obtained for a time-open , two-impulse transfer between fixed initial and 
final states. 


dJ = -H d 6 t 


(99) 


The classical condition for a time-open trajectory between fixed states is that the 
Hamiltonian function be equal to zero . Equation (99) indicates that the cost function 
J will be minimized when H 0. 


Orbit-to-Orbit Transfer 

Consider a transfer between two orbits, both fixed in shape and in orientation, 
for which the transfer time and the departure and arrival true anomalies 0 ^^ and 62 , 

respectively, are open (ref. 9) (fig. 10). What is the differential cost function for 
this solution? The differentials of the initial and final position and velocity vectors 
can be expressed as 


dR^ = dt^ 


dV^ = di^ 


dV2 = V 2 di2 


( 100 ) 


wherein t is a time measurement in the initial or final orbit . The differential of this 
time dr is related to the differential in true anomaly d 0 by the conservation of 
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Figure 10.- Transfer between two orbits. 


angular momentum (A.M.) . The magnitude of the A.M. vector in the two orbits may 
be written as 


h 


1 


h 


2 



r 

(101) i 


By solving for di^ and dt 2 from equation (101) and using the results in equation 
(100) , one obtains the position and velocity differentials 


dR 


1 




dV 


1 





( 102 ) 
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wherein n is the gravitational parameter and use has been made of the acceleration 
vector definition. Using equation (102) in equation (90) results in 


dJ 


h '>*1 " iT 




^2 




+ 






dt 


m 


(103) 


The coefficients of d0. and d0„ inequation (103) are 9J/90. and 9J/90„, 
respectively . ^ ^ 12 


Generalized Rendezvous 

As a final example , determine the differential cx>st function for the generalized 
rendezvous (G.R.) problem. A G.R. problem is characterized by an open-time 
transfer that is functionally dependent on the motion in the initial or final orbit. If 
the initial and final position and velocity vectors are given at the same time t^^ , then 

1. [Rj(tp, Vj(tj), t^] 

2. [R2(tp, V2(t^), tj] 

If they are also initiated with a transfer time 6t = t 2 ~ t^^, then the final position and 
velocity vectors are described as [R 2 (t 2 ) , V 2 (t 2 ) > t 2 l ■ If the transfer time is in- 
cremented by A6t to 6t' , then 6t' = 6t + A6t and the new final position and velocity 
vectors are described as [R2^^2 ^ * ^2 ^^2 ^ » ^2 ^ • Rewriting the gen- 

eral differential cost function (eq . (90) ) results in the following equation wherein 
the superscript ' refers to differentiation with respect to t. 


dJ = 


(xj^Vi - - 



+ 



(104) 
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In a G.R. problem, however, changes in time in the initial or final orbit t are equal 
to changes in the initial or final time , or 


= dtj 
*2 " '*‘2 


(105) 


Thus, dTj = Vj dtj et cetera . By rewriting equation (104) , making use of equa- 
tion (105), the following expression is obtained. 


dJ = - Vj-)dtj . - v/)dt2 

+ I X - A, pdR + (h ^ - H ]dt 
\m ml mlm m/m 

But 

(v/ - Vj-) = AV, 

(''2" - ''2') = '^''2 


and the impulses are in the direction of the multipliers 


AV^ = 


AVi 


AV2 = 


^^2 S 


(106) 


(107) 


(108) 


f 

It 




m 


r 

w 


Using equations (107) and (108) in equation (106) gives 


32 



(109) 
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as the differential cost function for the generalized rendezvous problem . Note that 
for optimal departure and arrival times from the initial and final orbits , the primer 
and its derivative must be orthogonal, or X 1 X (G .R . ) • This condition, as previously 
discussed , is exactly the condition on the primer vector derivative at an intermediate 
impulse point (eq. (32) for the nonimpulsive case) . A primer history for an optimal 
rendezvous with, for instance, three impulses may appear as shown in figure 11. 



Figure 11.- Optimal three-impulse 
rendezvous . 


That the primer and its derivative must also be orthogonal at an intermediate 
impulse point (t = t^) can be proved easily. On an optimal solution, the Hamiltonian 

function H must be continuous , or 


H =H 
m m 


T • - • -T - •T • 4- • +'p + 

X V - X V = X V - X V 
mm mm mm mm 


( 110 ) 


wherein X has been taken to be continuous . Now 
m 


V 

m 



( 111 ) 


Therefore, equation (110) resolves to 



( 112 ) 
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But , on an optimal solution , X is continuous , or 




T 


m 




= 0 


(113) 


But 


V 

m 


+ 


V 

m 


AV 


m 


m 


Thus , 


(114) 

s 


AV 


m 


X 


m 


m 


= 0 


(115) 


or X lx. 
m m 


INEQUALITY CONSTRAINTS IN PRIMER OPTIMAL N-IMPULSE SOLUTIONS 


The material in this section is based on reference 4. 


time t. 


Periapsis Constraint 

Consider the following sample problem . 

1. Depart from the initial position and velocity vectors and at the 
1‘ 

2. Arrive at the final position and velocity vectors R 2 and V 2 at the time t 2 - 


3. Problem: Find the optimum number of impulses and the minimum cost 
(eq. (47)) while requiring that the radius of periapsis of any transfer conic r of 

^i 

the N-impulse solution be larger than some value r , or r - r > 0 (fig. 12) . 

c Pj c 

To use the penalty function approach, consider a function 


Jp = pe 


1-p 


n 


(116) 


r 

k 


[ 
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Initial orbit 



Figure 12.- Radius-of-periapsis 
inequality constraint. 


wherein, as in figure 13, p is a positive constant, n is a suitably large positive 
integer, and p is 





(117) 


If the penalty cost for violating the constraint is designated as Jp and dJp/ dp 
as Jp' , then the following relationships apply . 

1 . When p » 1 , Jp — 0 and Jp' 2 ; 0 . 

2 . When p=l + e, Jp>0 and Jp' « 0 . 

3 . When 0<p<l, Jp>0 and Jp' < 0. 

The symbol e represents a small positive number , Note that Jp and Jp' remain 

finite for p = 1 because Jp and Jp' will be required to be continuous in the interval 

0 < p < 1 . Note that in this interval , Jp' < 0 , or that the parameter p will move in 

a direction that satisfies the constraint . If n is chosen sufficiently large , the pa- 
rameter p may be allowed to approach the constrained value as closely as desired . 
The parameter Jp is a single-walled function because p is constrained to take on 
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values only to the right of unity. A double-walled function (cost well) may be 
described by first creating a mirror function 


J 


q 


qe 


1-6 


m 


(118) 


F 

L 


wherein 


6 = a/n (119) 

s 

when a > 1 . If p = q and n = m , then the maximum width of the well is a - 1 as 
shown in figure 14. The cost well function is defined as 


Jx*, = J + J (120) 

W p q 




Figure 13.- Single-wall penalty 
function . 


Figure 14 . - Double-wall penalty 
function (cost well) . 


For a minimum penalty cost, 1 < p < a. 

The cost on an optimal N-impulse solution when inequality constraints are 
imposed on the transfer conics is 


J = J 


0 


+ J 


W 


( 121 ) 
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wherein 


and 


J 


0 


E 

^ — 1 


AV. 

1 


(122) 




N-1 

E 

i=l 




n 


Pie 


+ q.e 


1 - 6 . 

1 


m 


(123) 


For simplicity of nomenclature, consider a three-impulse solution (N = 3) . Then, 
the optimal differential cost can be written as 


->■’0 = " at ^ 


m 


m 




(124) 


wherein dSj and dS2 represent possible differential changes in the initial and 
final states (R,V,t). 

Because the penalty function is only a function of p, the differential is 


rtT rl 

^ ^ ajlT 


(125) 


wherein the subscripts 1 and 2 refer to the first and second conics , respectively . 
But p is a function of R and V only on the conics; because the constraint can be 
evaluated anywhere on these conics , it is convenient to choose the intermediate im- 
pulse time t . 

m 

Consider only the penalty cost on the first conic (dropping the subscripts 1 
and 2) . 


dJ,., = 


aj 


w dR + -51!- dV 


W 3p [3Rjj^ m 


av 


m 


m 


) 


(126) 
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The term dR requires no further change . The term dV must be expressed in 
m 

terms of the control parameters of the optimal N-impulse solution (dR^ ,dt^ ,dS^ , 
dS2> . Expanding dV^ results in 


av 

m 

~WR~ 

m 


dR + 
m 


(127) 


which may be written as 


dV = 6V + V dt 
m m mm 


(128) 


From the state transition matrix , the following expression may be obtained . 


/ X 

6V/tJ 


(129) 


Solving for results in 




6R, - (Pii 6R 

1 ^11 m 


(130) 


Now , the perturbations in the position vectors at t^^ and t^ are 


6R^^ = dRj - V^"^dt^ 


6R = dR - V dt 
m m mm 


Using equation (131) in equation (130) results in 


(131) 




dR. - v/dt, - 9ii(dR^ - dtj 
1 1 1 ll^m m m^ 


(132) 
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Using equation (132) in equation (128) results in the following equation for the differ 
ential dV 

m 

dV " = dR, - v/dt, - cp,,(dR - V "dt Ij 1 + V dt (133) 

m ^12 1 11 ^11\ m m mU mm 

The differential change in the penalty function (eq. (126)) for the first conic of the 
three-impulse solution is 


8J 


dJ.., = 


W 


w 


^ (p 

9R - ^12 

m av 

m 




dR 


m 


av 


m 


V + 
m 




(134) 


The total differential cost on an optimal three-impulse solution subject to an inequality 
constraint on the first conic is the sum of equations (124) and (134) . For various 
boundary conditions , the equations are modified as follows , 

1 . Initial departure time and position vector fixed 


dt^ = 0 
dR^ = 0 


(135) 


2 . Initial departure time fixed; initial true anomaly free 


dt 


1 


dR 


1 



(136) 


3. Generalized rendezvous 


dRj = dt^ 


(137) 
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The following is an example involving a radius-of-periapsis constraint . The 
magnitude of the radius of periapsis is 


r = h^/ (1 + e) (138) 

P 


wherein h and e are the magnitudes of the angular momentum H and the periapsis 
vector E , respectively . The parameter p is 





(139) L 


Rewriting equation (139) results in 


1 H^H 

^ 1 + |E| 


(140) 


wherein 


H = R X V (141) 

E = V X H - Pj^R/r (142) 

and is the gravitational parameter . Because H and E can be evaluated at 
t^ , the partial derivatives required in equation (134) can be expressed as 


-llL 

= ^ 

an 

+ ^ 

aE ’ 


3R 

m 

3H 

aR 

m 

aE 

aR 

m 

! (143) 

9P 

- ^ 

aH 

+ ^ 

aE 


av 

m 

an 

av 

m 

aE 

av 

m J 



r 
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r 
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Lighting Constraint 


As a final example , consider the problem of requiring an impulse to occur 
within a given lighting condition (fig. 15). The inertial frame (X,Y,Z) is located 
with respect to the Sun by the vector S . If the position vector R at the time of con- 
straint applicability is oriented with respect to the Sun vector by the angle \(/, then 
the function p for a lighting inequality constraint can be defined as u = \1//\1/ . , 

I- T/ 

Wherein and will describe the minimum and maximum bounds of the 

constraint. From figure 15, the angle vj/ may be expressed as 


= cos (144) 


wherein and fig are unit vectors in the R and S vector directions , 
respectively . 

The optimal impulsive solution requires not only the cost function J but also 
its gradient VJ. The total derivative of the cost function J is 


dJ 


= M ^ dR 
dp 3R 


(145) 


From the definitions given , the partial of p with respect to the vector R is 


dp _ ^S V R R / 

dR IRlvl/ . sin \1/ (146) 

• > ^mm ^ 


wherein I is the three-by-three identity matrix. From the conservation of angular 
momentum, the differential of R may be expressed as 


dR = 



(147) 


wherein 0 is the true anomaly and H is the angular momentum vector . Thus, the 
cost function J and its gradient can be readily expressed. The only question re- 
maining is how to determine the quadrant for the angle i|/ . 

If x)/ is measured positive by a counterclockwise rotation of the vector R into 
the vector S (fig. 15) , then the following test will define the quadrant for the 
angle i|/. 


1 i 1 < . 
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1. The expression v|/ = \[/, (R X V)'^(R X S) > 0 implies that the Sun is lead- 
ing the vector R . 

2. The expression \|> = 2x - \|/, (R X V)'^(R X S) < 0 implies that the Sun is 
lagging the vector R. 


r 

i 


Z 


s 



Figure 15.- Coordinate system for light- 
ing inequality constraint. 


CONCLUDING REMARKS 


r 

The necessary conditions that must be satisfied by the primer vector and its 
derivative on a general, two -body , N-impulse optimal trajectory have been developed. 
Example problems have been presented to illustrate a variety of equality and in- 
equality constraints on the solution . The Optimal , Multi , Impulse Rendezvous com- 
puter program was produced to encompass this general solution. 


Lyndon B . Johnson Space Center 

National Aeronautics and Space Administration 
Houston, Texas, May 15, 1975 
986-16-20-00-72 
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